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We consider a family of models describing the evolution under selection of a population whose 
dynamics can be related to the propagation of noisy traveling waves. For one particular model, that 
we shall call the exponential model, the properties of the traveling wave front can be calculated 
exactly, as well as the statistics of the genealogy of the population. One striking result is that, for 
this particular model, the genealogical trees have the same statistics as the trees of replicas in the 
' Parisi mean-field theory of spin glasses. We also find that in the exponential model, the coalescence 

, times along these trees grow like the logarithm of the population size. A phenomenological picture 

■ of the propagation of wave fronts that we introduced in a previous work, as well as our numerical 

| data, suggest that these statistics remain valid for a larger class of models, while the coalescence 

, | , times grow like the cube of the logarithm of the population size. 

Oh" 

^ ■ PACS numbers: 02.50.-r, 05.40.-a, 89. 75. He 

in : 

(N ; I. INTRODUCTION 

It has been recognized for a long time that there is a strong analogy between neo-darwinian evolution and statistical 
mechanics For an evolving population, there is an ongoing competition between the mutations which make 
J-j | individuals explore larger and larger regions of genome space and selection which tends to concentrate them at the 
* |*i i optimal fitness genomes. This is very similar to the competition between the energy and the entropy in statistical 
. ' mechanics. 

In the simplest models of evolution, one associates to each individual [1,0] (or to each species [J]) a single number 
which represents how fit this individual is to its environment. This fitness is transmitted to the offspring, up to small 
variations due to mutations. A higher fitness usually means a larger number of offspring [110,0,0,0,0, HI- If the size 
of the population is limited by the available resources, survivors are chosen at random among all the offspring. This 
leads in the long term to a selection effect: the descendants of individuals with low fitness are eliminated whereas the 
offspring of the individuals with high fitness tend to overrun the whole population. 

Our focus in this paper is a class of such models 0, 0, 0, 0, 0] describing the evolution of a population of fixed size 
N under asexual reproduction. The i-th individual is characterized by a single real number, Xi(g), which represents 
its adequacy to the environment. (This Xi(g) plays a role similar to fitness in the sense that offspring with higher 
Xi(g) will be selected; in the following, we shall simply call it the position of the individual.) At a generation g, the 
population is thus represented by a set of N real numbers Xi(g) for 1 < i < N. At each new generation, all individuals 
disappear and are replaced by some of their offspring: the j-th descendant of individual i has position Xi(g) + £i,j(g) 
. where £i,j(g) represents the effect of mutations from generation g to generation g + 1. Then comes the selection step: 
• ' at generation g + 1, one only keeps the N rightmost offspring among the descendants of all individuals at generation 
. g. One may consider two particular variants of this model: 

Model A: each individual has a fixed number k of offspring and all the ei,j(g) are independently distributed according 
to a given distribution p(e). For example, p(e) may be the uniform distribution between and 1. A realization of such 
r* ■ an evolution is shown in figure [T] Another example would be N branching random walks where the size of population 
is kept constant by eliminating the leftmost walk each time a branching event occurs. A visual representation of this 
latter example is shown in figure @. 

Model B: each individual has infinitely many offspring: the £i,j(g) are distributed according to a Poisson process of 
density tp(e) (this means that, with probabiliy tf>(e)de, there is one offspring of individual i with position between 
Xi{g) + e and Xi{g) + e + de). The density ip(e) is a priori arbitrary. The only constraints we impose are that ip{e) 
decays fast enough, when e increases, for the position not to diverge after one generation, and that /_ ip( e )de = oo, 
for the survival probability to be I. (This latter constraint implies in fact that each individual i has infinitely many 
offspring before the selection step.) 

As discussed in section [TlJ these models are related to noisy traveling wave equations ,of the Fisher-KPP type [lol . 
Til IT3 |. whic h ap pear in many contexts: disordered syste ms fl3L reaction-diffusion (l5l. [l6l. H?! fl8l| . fragmentation 
[lj| or QCD [la [Si, HI- A number of recent works [|, [H, [M0 US S H3, S HI, S3] focused on the fluctuations of 
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FIG. 1: Numerical simulation of the evolution of model A, with k — 2 and p(e) uniform between — ~ and \ for AT = 10. 
Upper plot: The filiation between each individual and its two offspring is shown. The individuals eliminated at each generation 
are shown by white circles, the surviving ones are shown by black disks. Lower plot: The noisy traveling wave front h g (x), 
constructed as in JTJ, is shown for the five generations of the upper plot. 




FIG. 2: A branching process for which the size N of the population is limited to 5. Each time the number of walks reaches 6, 
the leftmost walk is eliminated. Time goes downwards and the horizontal direction represents space. The actual population is 
represented in black, while the grey lines represent what the population would be for infinite N (i.e. in the absence of selection). 

the position these fronts, and this will allow us to predict how the fitness of the population evolves with the number 
of generations. 

Another interesting aspect of these models with stochastic evolution is their genealogy @ : one can associate to any 
group of individuals at a given generation its genealogical tree. One can then study how this tree fluctuates, and in 
particular what is the number of generations needed to reach their most recent common ancestor. The relationship 
between noisy traveling waves and genealogies is the main purpose of the present paper. 

While the models we consider here are difficult to solve for arbitrary p(e) and ip(e), one particular case of model B, 
with ^(e) = e _€ , turns out to be analytically solvable both for the statistics of the position of the population and for 
the properties of the genealogical trees. We shall call this case the "exponential model" and present its solution in 
section Ifffl 

As explained at the end of section [Til the exponential model is however non generic in the sense that it does not 
behave like a Fisher-KPP front. The generic case (which behaves like a noisy Fishcr-KPP equation but that we are 



not able to solve) and the exponential model can however be both described by a similar phcnomenological theory 
[1] , that we develop in section IIVI As a consequence, we argue that both the generic case and the exponential model 
have the same cumulants for the position of the front (up to a change of scale) , and that the genealogical trees have 
the same statistics in both models (up to a change of time scale). Numerical results, presented in section [Vl support 
these claims. 

II. THE LINK WITH NOISY FISHER-KPP FRONTS 

Our models arc nothing but stochastic models for the evolution of the positions of N individuals along the real 
axis. These positions form a cloud which does not spread: if an individual happens to fall far behind the cloud, it will 
have no surviving offspring, whereas the descendants of an individual far ahead of the cloud grow till they replace the 
whole population. With this picture in mind, it makes sense to describe the population by a front. Let Nh g (x) be 
the number of individuals with a position larger than x: 

I foo N 

h s( x )=M dz^8(z~ Xi(g)). (1) 

Jx »=i 

Clearly, h g (x) is a decreasing function with h g (—oo) = 1 and h g (+oo) = 0. In this section, we write the noisy equation 
which governs the evolution of this front. 

Let Nh* g+l (x) be the number of offspring on the right of x at generation g + 1 before the selection step. (So, for 
instance, h* +1 (—oo) is k in model A and oo in model B). Once h* +1 (x) is known, the selection step to get h g +i(x) is 
simply: 

h g +i(x) = min [l,h* g+1 (x)] . (2) 
Let us write the average and variance of h* +1 (x) for both models. 



A. Statistics of h* g+1 (x) for model A 



In model A, one can write 



JV 



Nh; +1 (x)=J2n g l) +1 (x), (3) 



where n g %i(x) is the total number of offspring before selection of the i-th individual of generation g which fall on the 
right of x. The probability that an offspring of i falls on the right of x is de p(e — a;,) and, as the k offspring of 
Xi(g) are independent, n^^x) has a binomial distribution. The average and variance are therefore given by 

n g+i( x ) = k J dep(e- Xi(g)), Variance (n^j^x)) = kj de p{e - Xi(g)) (l - j dep(e - Xi(g))j . (4) 

As the variables n^-yix) are uncorrelated, the average and variance of Nh* +1 (x) are simply from (|3|) the sums over 
i of the averages and variances of the n^+^x). For the average, one has 

/>OC />OC p 

Nh* g+1 (x) = kj deJ2p(e-x i (g))=-kJ de J dz p(e - z)Nh' g (z), (5) 
Where we used, from (fT|), 

N 

J2K X - X i(9))=-Nh' g (x). (6) 
Simplifying, and doing the same transformation for the variance, one finally gets 

/>oo 

1 - 2 / dz p(z) 



h* g+1 (x) = k I deh g (x - e)p(e), Variance (h* +1 (x)) = I deh g (x- e)p(e) 



for model A. (7) 

(Note that these average and variance are obtained for a given h g (x): they are not computed for the whole history. 



B. Statistics of h* +1 (x) for model B 



In model B, before the selection step, an individual at position Xi{g) has infinitely many offspring given by a 
Poisson process of density ip(x — Xi{gj). As Poisson processes are additive, the whole population (before selection) 
at generation g + 1 is also given by a Poisson process of density ^(x) with 

*(.t) = ip(x- Xl (g)) + --- + il>(x-x N (g)). (8) 
The number of individuals on the right of x is therefore a Poisson random number of average de ^(e), thus 

(9) 



Nh * g +i( x "> = Variance (Nh* g+1 (x)) 



deV(e) 



One can rewrite \I/(e) using the same trick as in ^ and ©. One finally gets 



h* g+1 (x) — J deh g (x — e)ij){e) and Variance (h* +1 (x)) ~ — j deh g (x — e)ip(e) for model B. (10) 



C. Front equations for both models and comparaison to Fisher-KPP fronts 
Comparing (fTUj) and {7}, one sees that one can write, for both models 



K+i( x ) = h *g+^ x ) + J h( x ) ^/Variance (h* g+1 (x)), 



where r] g (x) is a noise with r] g (x) = and Variance (r] g (x)^ = 1. Using ^ one finally gets 



h g+ i(x) = min 
h g +i(x) = min 



1, k / de h g (x — e)p(e) + 



1, / deh g (x - e)V(e) + 




k / deh g (x — e)p(e) 11 — 2/ cfz p(z) 



de h g (x — e)if)(e) 



(11) 

for model A, (12 A) 
for model B. (12 B) 



The precise distribution of r] g (x) depends on N and on the choice of the model. Far from both tips of the front, 
this distribution is Gaussian. At the tip, however, where h g (x) is of order 1/N, both h g (x) and its variance are 
comparable and the noise cannot be approximated by a Gaussian. (This is because the number of individuals is small 
and the discrete character of h g (x) cannot be forgotten anymore.) Furthermore, the noise is correlated in space but 
uncorrelated for different g. 

Thus, the precise expression of the noise r) g (x) is rather complicated, but its variance is 1, so that the amplitude of 
the whole noise term in (fT2")l decays as 1/yN as N becomes large. 
Equations (jTSJ) are very similar to the noisy Fisher-KPP equation: 



dh g (x) d 2 h g (x) 



dg 



dx 2 



h g (x) 



hg{ X f 



hg(x) 



hg(x) 2 , 



(13) 



where r]g(x) is a Gaussian noise with rj g {x) = and r]g(x)r]gi (x') = 8{g — g')S(x — x'). The noisy Fisher-KPP equation 
appears as a dual equation for the branching process A — > 2A (rate 1) and 2A — > A (rate 1/-/V) or, more simply, is an 
approximate equation valid for large N describing the fraction of A in the chemical reaction A + B — > 2A when the 
concentration of reactants is of order N [H, [3l|, [32| 

Comparing (fT2"|) and (| 13[) . the convolution of h g (x) by kp{e) or ip(e) in (fT2"|) spreads the front in the same way as the 
diffusion term in ()13j) . The same convolution induces the growth, similarly to the linear h g (x) term in (|13|) . as k<j)(e) 
and ip(e) both have an integral larger than 1. Thus, the fixed point h g (x) = is unstable. To balance the indefinite 
growth of h g (x), both (fT2"|) and (|13[) have a saturation mechanism (respectively the min(l, . . . ) and the —h g (x) 2 term) 
which makes h g (x) = 1 a stable fixed point. So, ignoring the noise terms (N — * oo), both (fT2"|) and Q13p describe a 
front which propagates from a stable phase h g (x) = 1 into an unstable phase h g (x) = 0. Finally, the noise terms in 
(fT2"|) and (fT3"|) have a similar amplitude of the order of y/h g (x) /N in the unstable region h g (x) <C 1. 

It is clear from the definitions of our models that the average velocity of the front is an increasing function of N. 
We first consider the limiting case N — > oo, which is equivalent to removing the noise term {rj g = 0) from (|12p and 



(|13p . To determine [12l ] the velocity of such traveling wave equations, it is usually sufficient to consider the linearized 
equation in the unstable region h g (x) <C 1 (where the saturation mechanism can be neglected). Looking for solutions 
of the form h g (x) ~ exp[— 7(2; — vg)], one gets a relation between the decay rate 7 and the velocity v = 11(7) that 
reads 



v(i) 



1 



In 



dep{e)e< t 



for model A, 



v(j) 



1 



In 



de^(e)e 7f 



for model B. 



(14) 



(For Fisher-KPP ([15]). one has ^(7) = 7 _1 + 7.) 

In many cases, when 11(7) is finite over some range of 7 and reaches a minimal value w(7o) for some finite positive 
decay rate 70, the selected velocity of the front for a steep enough initial condition [l2j is this minimal velocity w(7o)- 
For instance, for (fTS]) . one has 70 = 1 and the selected velocity is v(7o) = 2. Whenever this minimal velocity exists, 
we shall say that the model is in the universality class of the Fisher-KPP equation (fT5)) . For finite N, i.e. in presence 
of noise, there is a correction to this velocity and the front diffuses. We shall recall P in section IIVI that for the 
generic Fisher-KPP case, the correction to the velocity is of order 1 / In 2 N and that the diffusion constant is of order 
1/ In 3 N. 

There are however some choices of p{e) or tp( e ) f° r which v("f) is everywhere infinite or has no minimum. An 
example which wc study in some detail in section ITTTI is model B with ip(e) = e~ e , for which v(j) = 00 for all 7. We 
shall see among other things that, in presence of noise, the velocity of that front diverges as In In N for large N instead 
of converging to a finite value. 

It has been known for a long time that traveling wave equations are related to branching random walks [33l . |34| . 
This can be seen by considering a single individual at the origin at generation and by looking at the evolution of 
the probability Q g (x) that all of its descendants at generation g are on the left of x. In the case of model B with 
N = 00, one has 



Q g+ i(x) = ]J[1 - i/j(y)dy + if;(y)dyQ g (x - y)\ = exp (J 



dyip(y)(Q g (x — y) — I) 



(15) 



This equation describes the propagation of a front of the Fisher-KPP type, but where the unstable fixed point is 
at Q g = 1 instead of 0. For Q g close to 1, one gets exponentially decaying traveling wave solutions of the form 
1 — Qg( x ) exp[— 7(2; — vg)], with v = v(j) given by (|14|) for model B. (A similar calculation for model A leads to 
v(7) given by (fl4|) .) 

III. EXACT RESULTS FOR THE EXPONENTIAL MODEL 

In this section, we derive exact expressions (for large N) of the velocity, diffusion constant and coalescence times 
for model B with tp(e) = e~ £ . We first write some expressions valid for model B with an arbitrary density function 
V'(e), which we shall later apply to the exponential model. 

Before selection, the positions of the individuals at generation g + 1 are distributed according to a Poisson process 
of density *ff(x) defined in |(5J). We now wish to know the distribution of the A" rightmost individuals of this Poisson 
process, (i.e. of the offspring who survive the selection step.) We first consider the probability that there are no 
offspring on the right of x. Clearly, it is given by 



II [1 -*(*)<&] = exp (- / 



*(z) dz 



(16) 



Then, the probability that the rightmost offspring at generation g + 1 is in the interval [xi, x\ + dx\], and the second 
rightmost is in [x2 , X2 + dxi\ , up to the A^ + 1-st rightmost particle is 



^(xN+i)dxN+i ^(xN)dxN ■ ■ ■ ^(xi)dxi exp 



\l/(z) dz 



for xn+i < xn < ■■• < x\. 



(17) 



It will be more convenient not to specify the ordering of the A^ rightmost particles. Then the probability that the 
A^ + 1-st rightmost particle is in the interval [xn+i, xn+i + dx^+i] (as before) and that the A^ rightmost particles 
are in the intervals [xk, Xk + dxk] for 1 < k < N, with no constraint on the order of Xi, . . . , xjv, becomes 



1 



^(x N+ i)dx N+ i ty{x N )dx N 



^(xi)dxi exp — / 

y Jx N+1 



\J>(z) dz 



when xn+i < Xk for k = 1, . . . , N. (18) 



One obtains the probability that the TV+l-st rightmost particle is in the interval [xjv+i, XTv+i+cfccjv+i] by integrating 
over x±, . . . , xjj'- 



^(x N+ i)dx N+1 



poo 


N / oo \ 


/ $>(x)dx 


exp ( - / &(z) dz J 


J X N + 1 


\ Jx N+1 J 



(19) 



(As we imposed f^°° "0( e ) de = oo in the definition of the model, this distribution is normalized; see (8).) Finally, 
the probability of x%, . . . , xn given Xn+i is the ratio of ()18j) by (|19j) . One can see that, given the value of £jv+i) the 
distributions of x\{g + 1),.. . , xjvG? + 1) are independent and one gets that, given xjv+i, each of the N rightmost 
particles is in [x, x + dx] with probability 



f 



^{x)dx 
'°° *(x) dx 



for xn+i < x. 



(20) 



Therefore, to generate the whole population after selection at generation g + 1, one needs to calculate the density 
Hf(x) according to ([8]), then to choose the position of the N + 1-st rightmost particle according to (fT9|) and, finally, 
to generate independently the N rightmost particles xi(g + 1), . . . ,X]\[(g + 1) with the distribution (|20p . Note that 
the N + 1-st particle is not selected and is therefore eliminated after the N rightmost particles have been generated. 
This procedure is valid for any if>(e), but is in general complicated because © is not easy to handle analytically. 



A. Statistics of the position of the front in the exponential model 

In the exponential model ^(e) = e _e , however, everything becomes simpler: the Poisson process ((8|) becomes 

* exp (a;) = e -(.*-x a ) with X g = ln(^e Xl ^ +X2 ^ + - +XN ^, (21) 

which means that the offspring of the whole population is distributed as if they were the offspring of a single effective 
individual located at position X g . The distribution of the TV + 1-st rightmost particle (fT^|) becomes 

x N+1 =X g + z with Proba(z) = exp [-(N + l)z- e~ z ] , (22) 

and, once xn+i has been chosen, the distribution (f2"U|) of the Xk(g + 1) for k = 1, . . . , N becomes : 

x k (g + 1) = x N+1 + y k with Proba(y fc ) = eT Vk for y k > 0. (23) 

We now recall the calculation of the statistics of the position of the front [s[ which was done for a similar model 
in [IH, because we shall use later the same approach to calculate the statistics of the genealogical trees. 

There are many ways of defining the position of the front at a given generation g. One could consider the position 
of its center of mass, or the position of the rightmost or leftmost individual, or actually, any function of the positions 
Xk{g) such that a global shift of all the Xk{g) leads to the same shift in the position of the front. Because the front 
does not spread, the difference between two such definitions of the position does not grow with time so that, in the 
limit g — > oo, all these definitions lead to the same velocity, diffusion constant and higher cumulants. 

For the exponential model, it is convenient to use X g , defined in (|2ip . as the position of the front. Indeed, one can 
write 

AX g = X g+1 -X g = z + \n (e Vl + e V2 + ■ ■ ■ + e VN ) , (24) 

where the definitions and probability distributions of z and y k are given in (|22|) and |23|) . From (|24|) . the shifts AX g 
are uncorrclated random variables, and the average velocity vn and diffusion constant Djy of the front are given by 

v N = (AX g ), D N = (AX 2 g )-(AX g ) 2 . (25) 

More generally, all cumulants of the front position at a long time g are simply g times the cumulants of AX g . To 
compute theses cumulants, we evaluate the generating function G{(3) defined as 

e G(/3) = (e- pAX 3) = J dz Proba(z)e- /3z J dy x Proba(yi) • • • J dy N Proba(?Mr) {e Vl + ■ ■ ■ + e VN y , (26) 



and one obtains the cumulants by doing a small f3 expansion: 
Using ([2"2")l , the integral over z is easy: 



n>l 



J dz Proba(z)e~ 



A! 



dz exp + N + l)z — e~ z ] 



T(N + 1 + 0) 

r(jv + i) 



To calculate the integrals over yi in one can use the representation (valid for j3 > 0) 



r(/3) 7 



dX \ l3 ~ 1 e~ xz 



(27) 



(28) 



(29) 



with Z = e Vl + ■ ■ ■ + e yN . This leads to the factorization of the integrals over y%, - ■ ■ , J/jv- Replacing Proba(yfc) by its 
explicit expression from (|23| . one gets for j3 > (a similar calculation can be made for /?>—!) 



oGm r(N + i + (3) r~ dXX p-i L (X) N 



where 



f +oo 

J (A) = / dye' 



(30) 



(31) 



One can rewrite Iq (A) in several ways: 



A 



A 



W=A/ ^e-"-l + A(lnA + 7i j-l) + [e- A -(l-A)]-A / du 



1 - e 



n 



l + A(lnA + 7B-l)-^ 



(-1)* 



(32) 



fc+2 



fc=0 



(fc + l)(fc + 2)! 



where •je = _ r'(l) is the Euler constant. 

As Jo (A) is a monotonous decreasing function, the integral (j3"U|) is dominated by A close to 0. In fact, using (1301) , 
one can check that the range of values of A which dominate (|3l7|) is of the order of l/[Aln A]. Indeed, if one makes 
the change of variables 



ft = AAlnA, 



(33) 



one gets Jo (A) for values of fi of order 1: 

[I (\)] N ~exp[N\(hi\ + lE -l)), 



exp 



In A 



(hxfjt - In TV - In In A" + ~/ E - I] 



(34) 



e~M 1 + ft 



In ft - In In AT + -f E - 1 1 
ln~A + 2 



In /i — In In TV + 7s — 1 
rnA 



/'- 



where terms of order 1/A have been dropped. Replacing this expression into (|30p and using 

d/x^-^-^in^^—^r^), 



(35) 



one gets: 



1 



T{N + 1 + /?) 
T(iV + l)r(/3) (A In A)/ 3 
r(A + l + /3) 1 
r(AT + l) (A In A)' 3 



r(/J) 



r'(/3 + l) +T(/3+ 1)[- lnlnA + 7jB - 1] 



p fT'(p + l) 



in a V r(/3 + 1) 



In A 



- In In A + -f E - 1 



(36) 



(The next order is obtained in appendix ["/""J) The Stirling formula allows to simplify the expression: 

r(iV+1 + /?) 1 -1 + oflV (37) 



T(N+l) \N 
Then, one gets from (|36|) the following expression for the generating function 



gw = -em^N-JL ( lnln7V + ! _ 7£ _ + ( J_) . (38 ) 

(This expression was obtained assuming f3 > 0, but one can show that it remains valid for (3 > — 1 by using, instead 
of (|29|) . a different representation of .) Now one simply reads off the expressions of the cumulants of the position 
of the front by comparing the expansion of (|38p in powers of (3 and (|27l) : 



v N = = (AX g ) = lnlniV+ ^(lnln7V + l) + 



^.(fE..^ + (39) 

( X g) c = ( Ax n\ = n! CH = _n[_ y-J_ 

i>i 

up to terms of order lnlniV/ In 2 iV that are computed in appendix El The velocity vn diverges for large N, in contrast 
with models of the Fisher-KPP class for which wjv has a finite large N limit. Note that velocities which become infinite 
in the large N limit occur in other models of evolution with selection Q ■ 



B. Trees in the exponential model 



Let us now consider the ancestors of a group of p > 2 individuals chosen at random in the population (of size N). 
Looking at their genealogy, one observes a tree which fluctuates with the choice of the p individuals and which is 
characterized by its shape and coalescence times. 

For model B with an arbitrary density V"( e )j the probability of finding, at generation g + 1 before selection, an 
offspring in [x, x + dx] is ty(x) dx with ^ given by ((8]). On the other hand, the probability of finding in [x, x + dx] an 
offspring of Xi(g) is, by definition, ip(x — Xi(g)) dx. Therefore, given an offspring at generation g + 1 and position x, 
the probability that its parent was the i-th individual (at position Xi(g)) is 

W i{ x) = *(*-*;(*>) . (40) 
w(x) 

For general ip(e), these probabilities Wi(x) depend on x, making the calculation of these coalescence times difficult. 
In the exponential model, however, ([4*0)) becomes 

= e x i(9)-Xg = = (41) 

1 e xi(g) + . . . -{- e xjv(s) e s/i -| + e VN ' v ; 

where the — Xk{g) — XN + i(g) are the exponential variables of (|23p . Therefore the Wi do not depend on x. It 
follows that the probability q p that p individuals at generation g + 1 have the same ancestor at generation g is simply 

1p=(J2Wp\, (42) 

where the average is over the yi of (|4ip . After performing this average, all the terms in the sum over i become equal 
since the yt are identically distributed. Therefore 

r+ca /-+QO 

q p = N{Wf) =N dyx e- yi ■■■ dy N e - yN e PVl {e yi + ■ ■ ■ + e VN )~ p . (43) 



Using the representation (|29[) . one obtains 



N 



% = ^yyy j o dXX^I p (X)I (X) N -' (44) 
in terms of the function Iq(X) introduced in (|31[) and of its derivatives 

I p (X)=J dye^y~ XeV = (-)P-^-I (\) = \^p duuP~ 2 e- u . (45) 

For small A one has, to the leading order, 

/ (A)~l + A(lnA + 7B -l) ) / 1 (A)~-(lnA + 7B ) 1 I p (X) ~ for p > 2. (46) 

So far, (|44p is an exact expression and valid for arbitrary TV. From now on, we will work at leading order in In TV, 
leaving the extension to subleading orders to appendix El 

As for the obtention of ([55)) from the integral over A is dominated by the region where A is of order l/[jVln N]. 
Doing the same change of variable ^ = XNlnN, one gets Io(X) N — e _A1 and, using (|4*6"|) . X p ~ 1 I p (X) ~ (p — 2)!. 
Therefore, we obtain for p > 2 

Qp = 1 * 1 * ■ (47) 
In iV p — 1 

We see that for large N the probability that p branches merge is of the same order for all p, in contrast to the neutral 
model ([H, Hl| and appendix ICj) for which q p is of order 1/N V ~ 1 , so that q 2 ^S> (73 ^> (74 » • • • . 

To calculate the moments of the coalescence times, it is convenient to introduce the probability r p (k) that p randomly 
chosen individuals at generation g + 1 have exactly k ancestors at generation g. In one generation, at leading order in 
N, only a single coalescence may occur among the p individuals, and (|47[) tells us that the coalescence probability goes 
like 1/lniV (any additionnal coalescence at the same generation would in fact cost an additional power of 1/lnJV"; see 
appendix El) Consequently, we just need that p — k + 1 individuals coalesce to one ancestor, say individual number i 
(the probability is W[~ k+1 ), and that none of the other individuals have i as an ancestor (probability (1 — Wi)^ 1 ). 
Altogether, this reads 1 

r P (k) = ( fe ^) (j2wr k+1 (l . (48) 



The factor (1 — Wi) k 1 may be expanded and the average may be expressed with the help of the q p defined in 

r ^ k ) - (k P l) E (Y) H^-Vl- (49) 
Replacing (|4"T)) in (|4*5)l . one gets after some algebra 

rp{k) = hxN{p-k)(*p-k + l) ' (50) 

which holds for k < p. The probability r p (p) that there is no coalescence at all among the p individuals (that is to 
say, that all p have distinct ancestors) has a simple expression, which is obtained from a completeness relation: 

p-i _ , 

r p (p) = l-^r p (k)=l-?— . (51) 

fe=i 



1 In the mathematical litterature, one would rather use the transition rates Aj, q which give the probability that out of b individuals, the 
only event is the coalescence of the q first individuals [37I. l38l|. Clearly, r v {k) = L £^ A PjJ) _fe_|_i. All the Xf, q can be obtained through a 

measure A through A;, q = f} x q ~ 2 (l — x) b ~~ q A(dx). The exponential case corresponds to a uniform measure A, studied in |39| . 



The knowledge of the probabilities r p (k) in ((50]) and (f5Tj) allows one to determine (in the large TV limit) all the 
statistical properties of the trees. 

We introduce the probability P p (g) that p individuals have their first common ancestor a number of generations g 
in the past. For p > 2, one may write a recursion for P p {g) in the form 

p 

P P (g + 1) = XV^PfcG?) +r P (l)S° g . (52) 

k=2 

Using (|5T)|) and (f5Tj) . this becomes 

P p (g + 1) - P P (,) = -^P p ( 5 ) + g ^_£_ ft(ff ) + r , M . (53) 

In the large- N limit, the number of generations g over which the coalescence occurs is typically In N 3> 1 (since the 
coalescence probabilities scale like 1/lniV). It is then natural to introduce the rescaled variable t = g/lnN and the 
corresponding coalescence probability R p (t) dt = P p {g) dg. In this new variable, the recursion becomes for t > 

d -M> — (p.. w) t | > . t) ;., +1) w (54> 

This equation may be solved by introducing the generating function 

*(A,f)=^A p - 1 i? p (i), (55) 

which turns the summation over /c in (1541) into 



d 1 ^ d'fy 

— = [(1 - A) ln(l - A)]— - [ln(l - A)]*. (56) 



The general solution (which can be obtained by the method of characteristics) reads 



*(A,i) = T ^^(e- t ln(l-A)), (57) 



where </> is an arbitrary function. The initial condition for (|54p is the probability that all p individuals coalesce between 
times and dt (see (|3T)): 

R p (t = 0)dt = q p x^-dt = -^-j, (58) 

and thus, ([55]) becomes 

#(A,i = 0) = -ln(l-A). (59) 

This leads to 

*(A,t) = |(l-A) e_t - 1 . (60) 

The expansion of (|6"0|) in powers of A using 



leads through ([55)1 to 



r(a)^jr(p + i)- 



which is just a polynomial of order p — 1 in the variable e *. More explicitely, for the first values of p, one finds 



R 2 (t)=e~ t , R 3 (t) = -e ~* - e 
The average coalescence times (using ([62]) ) arc 

00 Z'+oo poo 

(T P ) = V ffP P (<?) =lnN dttRp{t)=lnN dt 

n-n JO JO 



\ i24(t) = ye-* 



2e~ 



1 

2 6 



-3t 



3=0 

and one gets 



l-(l-e-*) 1- 



1 - 



1 - 



p-1 



(T 2 )=lniV ) (T 3 ) = -(T 2 ), (T 4 ) = -(T 2 >, 



(63) 



(64) 



(65) 



These expressions contrast with a neutral model of coalescence with no selection [3a . l40l ] where at each generation 
one would choose the N survivors at random among all the offspring at generation g + 1 (see appendix [Cj) : 



(T 2 ncutral ) =0(N) , (T^ lcutral ) = l(T 2 noutral ) , (T 4 

3 



neutral 1 



ncutralx 



(66) 



(Table [I] compares the frequencies of the trees in the cases with and without selection. 





Neutral case 


Exponential model 
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Neutral case 



Exponential model 



TABLE I: Probabilities of observing each of the possible genealogical trees for three and four individuals in the neutral case 
and in the exponential model 



As shown in appendix [Bj the ratios ([rJS"|) arc on the other hand identical to those which would be computed if the 
genealogical trees had the same statistical properties as mean-field spin glasses [3^, |4l| . 

We also see that (T p ) in (f6"5| scales like In AT for any fixed value of p, which means that on average, a given number 
of individuals have their first common ancestor at of order lnA^ generations in the past. It is however interesting to 
note that for large p, 

(Tp)~lnATxmlnp (67) 

which is obtained by using, from ([B"2"|) . R p {t) ~ jp" cxp '"'' ~ rf e -exp[-(t-lninp)] f or i ar g e p ; R p (t) becomes a Gumbel 
distribution of width of order 1 centered at In In p. 



IV. PHENOMENOLOGICAL EXTENSION TO GENERIC MODELS 



The exponential model had the advantage of being exactly solvable, but as already mentioned, it is non-generic 
because the velocity — > 00 as N — > 00, in contrast to models of the Fisher-KPP type. We do not know how to 
calculate directly the velocity vn, diffusion constant or the coalescence times of the generic Fisher-KPP case. 
One can however use a phenomenological picture of front propagation [8[ and ancestry, which is fully consistent with 
exact calculations in the case of the exponential model, and with numerical simulations in the generic case. 



A. Picture of the propagation of fluctuating pulled fronts 



Let us recall the phenomenological picture of front propagation which emerged from [g, |42|]. In this picture, most 
of the time, the front evolves in a deterministic way well reproduced by an equation obtained from (|12|) by removing 
the noise term, and by adding a cutoff which takes into account the discreteness of the number of individuals: This 
ensures that h g (x) cannot take values less than 1/N. The evolution equation in the case of model B reads [42| 



, min 1, / deib(e)h n (x — e) I if that number is larger than 1/N 
hg+i(x) = ^ V J 7 _ (68) 

^ otherwise. 

In the exponential model {ip(e) = e~ e ), it is easy to see that the solution to ((G8l) is 



(69) 




where the parameter Y g can be used as the definition of the position of the front. Substituting ([B"9"|) into 
obtains the velocity 



"cutoff 



Y g+1 -Y g = ln(ln N + 1) ~ In In N, (70) 



which does agree, to leading order, with the exact expression ([39]) . 

For a front in the Fishcr-KPP class, the cutoff theory can also be worked out [42|. One obtains 

h g (x) « L sin (*^) ^° {X - Ya) ^d = Y g+1 — Y g ~ v( 7o ) (71) 

where 17(7) is given by (|14p . 70 is the value of 7 which minimizes ^(7), and Lq = (]nN)/yo is the length of the front, 
from the region where h g is of order 1 to the region where it cancels. The expression of h g (x) in (|7ip is only valid for 
h g (x) -C 1 and x — Y g < Lq. 

By convention, we shall define 70 = 1 in the exponential case. Then, both in (|69|) and in ([7T1) . the front has 
essentially an exponential decay with rate 70 and its length is Lq = (\nN)/jo. 

So far, (|70p and (|7ip have been obtained from a purely deterministic calculation, where only the discreteness of 
h g (x) has been taken into account. Stochasticity may be put back in the picture for the generic (Fishcr-KPP) case 
in the following way [1] : 

From time to time, a rare fluctuation sends a few individuals ahead of the front at a distance 5 from its tip. This 
occurs during the time interval dt with a probability p(S) d5 dt where p(5) was assumed Q to be 

p(S) = C ie -^ 5 (72) 

for <5 large enough. C\ is a given constant. 

These individuals then multiply and build up their own front in an essentially deterministic way. After about Lq 
generations, the descendants of these individuals have mixed up with the individuals that stem from the rest of the 
front. The effect of this rare fluctuation is therefore to pull ahead the front by a quantity R(5) which, in the generic 
(Fisher-KPP) case, is given Q by 

70 S " 

(73) 



1 / e 7 ° 5 \ 



where C2 is another constant and a = 3. Finally, in Q it was argued that 

C X C 2 ^ 2 7oW "( 7o ). (74) 

As we shall show in the next section, the same picture applies to the exponential model with some slight modifica- 
tions: in (|73[) . one needs to take a = 1 instead of a — 3, everywhere 70 must be replaced by 1, one should replace Q74|) 
by C\ = C2 = 1 and the relaxation time of a fluctuation by 1 instead of Lq. 

With these ingredients, it is not difficult to write the generating function of the position Y g of the front: 

(e"^) - e 9G ^ where G(J3) = -(3v cut oS + / dS p{6) (e^^ - l) . (75) 



The first term in G(/3) is due to the deterministic motion, while the integral represents the effect of the forward rare 
fluctuations. In the case of the exponential model, this expression leads to (j3l?l) . up to terms of order 1/lniV for 
the velocity and of order In In N/ In 2 A" for the other cumulants. In the generic Fisher-KPP case, the average front 
velocity, diffusion constant and higher order cumulants are found from ([75)) to be Q 

, v ^ 2 7 V'(7o) , a/// \ 2 3 In In ^7 2 «"(7o) 
v N = w(7o) -z 2 AT 1" lo v (7o)tt : g Ar H = "(To) - 



2 In 2 TV '° w ' ; In 3 TV Wl " 2(lnJV + 3 In In TV) 2 

Dn = 7ov"(7o)-K- + ---, (76) 
3 In N 

<? In TV 

One important aspect of (|73[) is that when <5 is of order (alnLo)/70j the front is shifted by one additional unit in 
position due to this fluctuation. This means that a large fraction of the population is replaced by the descendants of 
the individuals produced by this fluctuation. Thus, when one considers a given number of individuals at generation 
<7, the most probable is that their most recent common ancestor belongs to one of these fluctuations that triggered 
shifts of order 1 in the position of the front in the past generations. According to (|72p . such events occur once every 
Ag ~ Lq generations. Ag is likely to give the order of magnitude of the average coalescence times. In scction llV C] we 
shall build on this to obtain the statistics of the genealogical trees and the coalescence times in the generic Fisher-KPP 
case. But first, we show that this phenomenological picture is consistent with the exact results ([55)1 for the exponential 
model. 



B. Exponential case 

Since the exponential model can be solved exactly (section IIII|) , one can test in this case our phenomenological 
picture of section IIVAI Let us first show that (f72|) gives the correct distribution of fluctuations. 

In the exponential case at any generation g, the front is built according to (|23[) by drawing A~ independent expo- 
nential random numbers yk which represent the positions of the particles relative to a common origin xjv+i- There is 
a probability (1 — e~ v ) N that none of the yk are on the right of y; therefore the distribution of the rightmost yk is 

Proba(y rightmost ) = N (l - e^w^c^" 1 e -xw**»™. ~ e xp [- (y right most - hiN) - e -(fcw*»-*-inAr)j _ (77) 

^rightmost is the distance between the rightmost particle and the N + 1-st rightmost particle (before selection). We 
define the length I of the front as I = ^rightmost ■ (A more natural definition could have been the distance between the 
rightmost and the leftmost particles, which is obtained by replacing A^ by A^ — 1 in the previous equation. For large N, 
this difference between these two definitions is negligible.) The average length of the front is therefore (I) ~ In N + -fE 
with fluctuations of order 1 given by a Gumbel distribution, and the probability to observe a large fluctuation where 
I = \nN + S with S 3> 1 is given by 

p(5) ~ exp [-S - e~ s ] ~ exp [-S] , (78) 

which is the same as (|72|) . 

We now wish to know the effect of such a fluctuation on the position of the front. As the shape of the front is 
decorrelated between two successive generations, the relaxation time of a fluctuation is 1 and it is sufficient to compute 
AX g given the value of 5 at generation g. Given the value of I = ^rightmost, the distribution (|2"3"|) of the A^ — 1 other 
yk become 

e ~Vk 

Proba(y fc ) = j for < y k < I. (79) 

1 — e 1 

As in (|2l)|) . we introduce the generating function of the displacement AX g given the value of I: 



(e" MX »|l) = j dz Proba(z)e-^ / d Vl Proba( yi ) ■■■ / dy N ^ Proba(yiv-i) {e Vl + ■■■ + e^- 1 + e l ) , 

I" d Vl e-^ ■ ■ ■ [ dy N - X e"^- 1 (e^ +■■■ + e^" 1 + eJV' 3 , 



r(iv + i) (i - e-^- 1 Jo 



(80) 







where ((2"5} and ([7U|) were used. By using the same representation PS]) that led to (|30p . one gets 



r(iv + i + /3) 



r(jv + i)rcs) y 

which, in terms of Iq(X) defined in (|31[) . is the same as 

T(N + 1 + /3) 



1-e" 



-y-Ae« 



iV-1 



-Ae' 



(81) 



r(Jv + iTO) y 



dW 13 - 



I (X)-e- l I Q (Xe l ) 



1 - e" 



AT-l 



-Ae' 



where, using (f32l) . 



/ (A) - e-'/oCAe') 



1 -e" 



1 — A- 



1 - e" 



+ 00 

E 



(-1)* 



(k + l)(fc + 2)! 



fc+2 



J(fc+1) _ 1 

1 - e~ l 



(82) 



(83) 



Expressions (fB"2")) and are valid for any value of I . We now consider a large fluctuation Z = In N + 8 with 
1 <C 8 < lnlnAf. As for (|30[) . the integral is dominated by values of A of order l/[7VlniV]. Making as before the 
change of variable ft = \N\nN, and dropping all the terms of order 1/N, one gets 



/ (A)-e- ; / (A e ') 
l-e~ l 



N-l 



exp 



-fill 



In TV 



E 

fc=0 



(-1)* 



(k+ l)(fc + 2)! VlniV 



fc+2 



,<S(fc+l) 



(84) 



We are only interested in the leading order in 1/lniV. Dropping higher order terms, one gets, in 



-PAX g 



S) 



T(N + l+f3) 1_ 

T(N+1)T(J3) [NliiN}? J 



dfi jjL^ 1 exp 



8 + e 6 
InN 



[In NY 



1 



In TV 



(85) 



where ([57} has been used and where 8 was neglected compared to e 5 . 

This means that up to the order l/[lniV] we are considering, AX g given 8 is deterministic with 



f cutoff + R(8), 



AX g (8) ~ In In AT + In 1 



In TV 



(86) 



where we used ([701) and (J73J) with C2 = a = 70 = 1. 

The phenomenological picture we developed for the generic case is therefore justified for the exponential case: each 
rare fluctuation of size 8 in the length of the front leads to a shift R(8), given by (|73|) . for the position of the front. 



C. Genealogical trees 

With the above scenario, one can also build a simplified picture for the evolution of a population. We assume that, 
at each generation, there is with a small probability a fluctuation of amplitude / produced by an individual ahead of 
the front. The long term effect of this fluctuation is that a fraction / of the population is replaced by the descendants 
of this individual. 

One can now relate the probability distribution of / to the phenomenological picture of front propagation. Starting 
with a front at position Y go at generation go, we consider its position Y g at a generation g > go- If no important 
fluctuation has occurred, the tail of the front is given by 

hno fluctuation^, g) oc ( X ~ Y ^ """") with Yg° = Y go + v cutoS (g - go). (87) 

(See (|7ip: for simplicity, we neglect the Sine prefactor in the tail as it is a slow varying factor which, to the leading 
order, does not change our final result.) 

If instead a fluctuation has occurred, generated by an individual ahead of the front by a distance 8, then the shape 
is eventually described by 

/^fluctuation^, <?) CX e^To (x-yj— ) ^ ^fluctuation = y^ + _ g ^ + R ( § y ( g g) 



FIG. 3: Effect of a fluctuation of a front. The dashed line is the front (|87[) in the absence of a fluctuation. The plain line is 
the front ()88|) if a rare fluctuation occured. The grey area represents the contribution to the front from the descendants of the 
fluctuation. After the front has relaxed, they represent a proportion / of the whole population. 



that is, the front is pulled ahead by R(S). If one assumes that the extra mass in the front with fluctuation (in grey 
in figure [3J is due to the fraction / of descendants originating from the fluctuation, then one gets h no fluctuation = 
(1 — /^fluctuation- The substitution of (|57|) and ([55)1 yields 



- 1 _ p -7ofl(<5) 



(89) 



This equation defines the mapping between the / and the S representations of the phenomenological model. The 
probability distribution of S in (j?2")) and the expression ([75)) of R(S) implies the following distribution of /: 



Proba(/) = 



CiC 2 1 



(90) 



(Note that this expression cannot be valid down to / = for the distribution to be normalized. One should therefore 
consider that ([9"0"1) is valid above a certain small threshold /mm- This threshold has no effect on the correlations 
calculated below.) 

Using (|7^|) and a = 3 in the Fisher-KPP case, and C\ = C2 = 70 = ce — 1 in the exponential case (see section HVB I) , 
one gets 



Proba(/) 



1 1 

lnJV/2 

^ 2 7oV'(7o) J_ 
P 



In 3 N 



for the exponential model, 

for the generic Fisher-KPP case. 



(91) 



In this model, p individuals may coalesce if they belong to the fraction / of individuals that are the descendants of a 
fluctuation. The probability of such an event thus reads 



q p = f d/Proba(/)/^ 
Jo 



dC 2 1 
70^0 P - 1 



(92) 



which, for the exponential model, is identical to the exact asymptotic result in (|47|) . 

The coalescence probabilities in one generation r p (k) may be obtained in a straightforward way in this model. One 
first chooses the k — 1 individuals among p that do not have a common ancestor in the previous generation. The 
latter have to be part of the fraction 1 — / of individuals, while the remaining p — k + 1 individuals that have their 
common ancestor in the previous generation must belong to the fraction /. Thus 



r p (k) 



P 

k - 1 



df Proba(/)/ 



p-fc+i 



(1 - if- 1 = 



p 



7o l« { P -k)(j)-k + iy 



(93) 



with the same result as in (|50|) for the exponential model. 2 At this point, the combinatorics to get the coalescence 
probabilities and average times are the same as in the exact calculation for the exponential model in section ["HI Bl So, 
for the exponential model we recover the results of section IIII Bl and for the generic Fishcr-KPP case, we get instead 

1" 7o u (7o) 

while the ratios (Tj) / '(T2) are the same (]65p as for the exponential model, in agreement with the results of numerical 
simulations of @ and of section fVl below. Indeed, the r p (fc)'s given in ([50)1 and (|93|) are identical except for an overall 
constant which cancels out in the ratios. 

We note an interesting relation between the average coalescence time and the front diffusion constant, valid both 
in the exponential model and in the generic Fishcr-KPP case: 

D N x (T 2 ) ~ * (95) 
J 7o 

We will test numerically this identity in section fVl 

As a side remark, we note that if Proba(J) of (f9"Tj) is replaced by Cste/~ a with a — > 3 (instead of a = 2 in our 
selective evolution models), then the ratios of the coalescence times are identical to those obtained for evolution 
models without selection, see appendix [C] 



V. NUMERICAL SIMULATIONS 



A. Algorithms 



In order to measure the velocity and diffusion constant of our models, it is sufficient to follow the evolution of the 
positions of the individuals. In the case of model A, at each generation, one first draws at random the k offspring of 
each individual and then one keeps the N rightmost offspring as the new population. This can be done in a computer 
time linear in N. For model B, one can start by drawing at random the two rightmost offspring of each individual. 
If Z is the position of the N-th. rightmost offspring out of this first set of 2N, then one draws for each individual all 
its remaining offspring which are larger than Z. Then, taking the N rightmost individuals among those drawn gives 
the new population. 

We measured the diffusion constants Dn as in [43[, using Dn = ((X go+g — X go — vg) 2 )/g for a large g. (This 
expression is in principle only valid in the g — * 00 limit.) In practice, we have to choose an appropriate value of g and 
average over many runs. For each value of N, we measured the diffusion constant twice, once with g ~ 2 In 3 N and 
once with g ps 10 In 3 N, and we have plotted both values with the same symbol. The fact that one cannot distinguish 
the two sets of data indicates that the values of g we took are large enough and that we accumulated enough statistics. 

To measure the statistics of the genealogical trees in the population, one needs to memorize more information than 
simply the positions of the individuals in the current generation. The most naive method would be to record the 
whole history of the population, keeping for all individuals in all generations their positions and parents, and then to 
analyze at the end the whole genealogical tree. This is clearly too time and memory consuming. Instead, we used the 
three following algorithms. 

The first algorithm consists in working with a matrix T g , the element T g (i,j) being the age of the most recent 
common ancestor of the pair of individuals i and j at generation g. This matrix is simple to update: if j and j' are 
the parents of i and i', then T g+ \(i 1 i') = 1+T g (j, j') for i ^ i' and T g+ i(i,i) = 0. By sampling random elements of the 
matrix at different generations, one obtains the average value of the coalescence time between two individuals. The 
nice thing is that, due to the ultrametric structure of the tree (for any i, j and k, T g (i, j) < max [T g (i, k), T g (j, ), 
no more information is needed to compute the coalescence times of three or more individuals: the age of the most 
recent ancestor of p individuals i\, . . . , i p is simply given by max \T g (ii, T g (ii, 13), . . . , T g (i\, . This method 
is appropriate for values of N up to about 10 3 as it takes a long time of order A^ 2 to update the matrix at each 
generation. 

In the second algorithm, instead of working with this matrix T g (i,j), we take advantage of the tree structure of 
the genealogy by recording only its "relevant" nodes: at generation g, we say that a node is "relevant" if it is an 



2 In the language of the transition rates A M defined in [HIH, one would write X b q = df p(f)f{l - f) b ~ q oc /* df f q ~ 2 {l - f) b ~ q . 
It is the A-coalescent with the uniform measure, i.e. the Bolthausen-Sznitman coalescent. 



individual of the current generation g or if it is the first common ancestor of any pair of individuals of the current 
generation. Clearly, the "relevant" nodes have a tree structure (the first common ancestor of any two "relevant" nodes 
is a "relevant" node), which we record as well. The leaves of this tree are the current generation, and the root is the 
most recent common ancestor of the whole population. This tree is simple to update: if, after one timestep, a node 
has no child, it is removed and its parent is updated. If a node has only one child, it is removed as well and its child 
and parent get directly connected. If the root of the tree has only one child, it is removed and its child becomes the 
new root. As can be seen easily, the tree has at most 2N — 1 nodes and it can be updated in a time of order N. The 
extraction of the interesting information from the tree is also very fast: if a node has p children, and these children 
are the ancestors of a\, . . . , a p individuals of the current generation, then this node is the most common ancestor of 
J2iytj OLiOLj pairs of individuals. More generally, this node is the most common ancestor of (^~' i q ai ) — J2i groups of 
q individuals in the current generation. By computing this quantity on each node of the tree, one obtains the average 
(or even the distribution) of all the coalescence times within the current generation in a computer time of order N. 
This algorithm turns out to be very fast and we used it for N up to about 10 6 . 

The third algorithm only works for a limited class of models, for which the positions Xi(g) are integers: instead of 
recording the N positions, one only needs to record the number of individuals at a given site. The typical width of 
the front and, therefore, the number of variables to handle, are only of order ]nN. Let us, then, consider model B 
with ^>(e) given as a sum of Dirac functions: ip{e) = Ylq 4>q5{t — <?)• This means that, before selection, an individual 
at position x has a number of offspring at position x + q which has a Poisson distribution of average <f> q . Considering 
now the whole population, the number of offspring at time g + 1 and site y is also a random Poisson number of average 
^2 x n(x,g)(j)y- x , where n{x,g) is the number of individuals at site x and time g (compare to ©)■ To simplify, we 
consider only cases where <j> q = for q larger than some qo, so that one can easily update the system from right to 
left by drawing Poisson numbers and stopping when the total number of individuals at time g + 1 reaches N. So far, 
the method described allows us to update the positions of the particles, and therefore to extract the velocity and the 
diffusion constant, in a time proportional to In TV per generation. A similar method has already been used in [42l l43j 
to simulate populations up to N ~ 10 100 . To extract the coalescence times, one needs to keep more information. 
The difficulty resides in the fact that the many individuals at a given position usually have different ancestors. 
To overcome this difficulty, one can consider the average coalescence times T g (x,x / ) of two different individuals at 
respective positions x and x' . To update that matrix, one starts from the probability that an individual of generation 
g + 1 and position y is the offspring of an individual who was at position x: 

Proba(y comes from x) = — — ^li^l x _ (95) 

Ex' n O ',9)<t>y-x> 

(Compare to (00]).) Then, one obtains that 

/ 5*' \ 

T g+ i(y, y') = 1 + Proba(y comes from x) Proba(y' comes from x')T g (x, x 1 ) I 1 v- 2 — r ) . (97) 

(The term in parenthesis is the probability that individuals at positions y and y' come from two different parents 
given the parents' positions x and x'.) Then, the average coalescence time of two individuals in the population is 
simply given by 

Therefore, by storing a matrix of size In 2 N which can be updated in a time In 4 N, one can obtain the average 
coalescence time of two individuals. An interesting observation is that this algorithm simulates one possible realization 
of the positions of the particles; however, the quantity T g (x, y) is actually an average over all the possible genealogical 
trees in the population given that realization of the positions over time of the particles. A complexity in time of order 
In 4 N allows already to simulate rather large systems. However, a further optimization is possible in the special case 
where <f> q is constant for q < qg. For that specific model, additional simplifications occur (one can write a recursion 
on the matrix elements) and the matrix T g (x, x') can be updated in a time of only In 2 N. This allows one to study 
systems of size N up to about I0 50 in a few weeks time on standard desktop computers. There is, unfortunately, 
not enough information in the matrix T g (x, x') to extract the average coalescence time of three (or more) individuals: 
to that purpose, one needs to simulate a tensor with three (or more) indices which can be updated with rules very 
similar to (|97p . Because of this extra complexity, we only measured the average coalescence time of three individuals 
for values of N up to I0 20 . 



B. Results 



Using this last algorithm, we have simulated model B for ip(e) = i Yl n <o — n ) U P to AT = 10 50 . The velocity 
and diffusion constants are shown in figure [H compared to the predictions ((75|) in plain lines. There is still a small 
visible difference between numerics and theory, but this difference gets smaller as N increases. In order to obtain a 
better fit, we have included sublcading corrections by changing the denominator (In AT + 31nlniV) 2 for the velocity 
in ([75]) into (In ./V + 3 In In TV — 3.5) 2 . Similarly, we changed the denominator (lnA^) 3 for the diffusion constant in ([To]) 
into (In A" + 3 In In TV — 3.5) 3 . With these subleading terms (in dotted lines on the figure), the fit is almost perfect 
over more than 40 orders of magnitude. 




(In N + 3 In In N - 3.5) 2 (In N + 3 In In N - 3.5) 3 



10 3 10 5 10 10 10 30 10 50 

FIG. 4: Numerical simulations of model B with ij){e) = \ ^2 n<0 S(e — n). The circles are the correction to the velocity and 
the triangles the diffusion constant, as a function of N. The plain lines are the predictions (|76p . The dotted lines are the 
predictions (|76[1 with, for both quantities, the same subleading terms added in the denominators. (The scale on the N axis is 
proportional to lnlniV.) 

We have no theory to justify these extra subleading terms, but we simply notice that it is possible to fit both the 
correction to the velocity and the diffusion constant using the same subleading terms in the denominators of their 
respective expressions. 

For the same model, (T?) is shown on figure [5] (using circles), compared to the prediction (|94p in plain lines. As 
for the velocity and diffusion constant, there is still a small visible difference and we obtain a better fit if we include 
subleading terms (in dotted lines): guided by ([55)1 and the fit used for the diffusion constant in figured! we changed 
the numerator of (JMJ) from (In N) 3 into (In N + 3 In In AT - 3.5) 3 . On the same figure, (T2) for the exponential model 
is shown (using triangles), compared with the exact prediction (|65p (T2) ~ In AT. Here again, the fit is improved by 
including the subleading corrections (|A15j) (T2) ~ In A 7 " + In In A" obtained in appendix [Al 

Figure [5] combines data from figures 0] and [5] The triangles are the ratio of the diffusion constant and of the 
correction to the velocity to the power 3/2. For large Af, this should converge to a constant which we can compute 
from (|76j). The circles are the product of the diffusion constant and of the coalescence time (T2), which we expect to 
converge to the value given in (|95p . The horizontal lines on the figure represent both predictions. 

Finally, figure [7] shows the ratio (T^)/{T2) as a function of A" up to TV = 10 20 . The ratio is very close to 1.25 for 
large AT, which is the prediction of the phenomenological theory of section HV CI (see also (jM}). 

VI. CONCLUSION 

In the present work, we have solved exactly a simple model of evolution with selection, the exponential model 
of section IIIII For this model, we have calculated the velocity and the diffusion constant (|39]l of the parameter 
representing the adequacy of the population to its environment, as well as the coalescence times which characterize 
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FIG. 5: Numerical simulations of (T2) for model B with ij){e) — jE n <o''( e — n ) (circles) and for the exponential model 
(triangles). The plain lines are the predictions (|65p for (T2) and (|94[) . while the dotted lines are the same predictions with 
some subleading term: for the generic case, we used subleading terms suggested by (|95[1 and the fit of figure [11 and for the 
exponential model the exact results (|A15|I . 
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FIG. 6: Numerical simulations of model B with tp(e) = t Yl n <o 8{e — n). The circles represent the product Dn x (T2) compared 
to the prediction (|95p . The triangles are the ratio of the diffusion constant and the correction to the velocity to the power 3/2, 
compared to tt^8/v" (70 )/(3yq), which is the prediction obtained from ([76jl . 



the genealogy. We have shown that the statistical properties of the genealogical trees are identical to those trees 
which appear in the Parisi mean field theory of spin glasses [HI, [4{| ■ They therefore follow the Bolthausen Sznitman 
statistics [H, EH , in contrast to the case of evolution without selection which obeys the statistics of the Kingman 
coalcscent. 

The reason why the exponential model is exactly soluble is that, going from one generation to the next, the only 
relevant information on the position of the individuals is contained in one single variable X g defined in (|21j) . The 
exponential model belongs to a larger class of models parametrized by a single function p (for model A) or -0 (for 
model B). We have not been able to solve the generic case and, unfortunately, the exponential model is special: 
while the generic case can be described by a Fisher-KPP front, with a velocity which converges when N — > 00, the 
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FIG. 7: Numerical simulations of model B with ip( e ) = i X^n<o — n )- The circles represent the ratio (T3) / '{T2) as a function 
of N, compared to the result 5/4 suggested by the phenomenological theory of section TlV CI (The scale on the N axis is 
proportional to 1/lniV.) 

velocity of the front associated to the exponential model diverges when N — > 00. We have however constructed a 
phenomenological picture of front propagation which can be used both for the exponential model and for the generic 
Fisher-KPP case, and which also provides predictions for the genealogy. Within this picture, we have that the average 
coalescence times scale like In 3 N with the size TV of the population for the generic Fisher-KPP case (while it grows 
like In N for the exponential model) , and that the structure of the trees is the same as in the Parisi mean- field theory 
of spin glasses. 

Proving the validity of the phenomenological picture for generic models is an interesting open question for future 
research. Understanding more deeply why our models of selective evolution are related to spin glasses would also 
deserve some efforts. Lastly, it would be interesting to study genealogies in other models of selective evolution to 
test the robustness of our results. 
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APPENDIX A: EXACT RESULTS FOR THE EXPONENTIAL MODEL INCLUDING SUBLEADING 

ORDERS 

In this appendix, we obtain higher orders in the large IniV expansion, for the statistics of the position of the front 
and for the coalescence probabilities in the exponential model. 

1. Front position statistics 

The exact expression for the cumulants of the front velocity was given in (|30|) in terms of the function Iq defined 
in (j3l"j) . Discarding all the terms of order 1/N or smaller, one can use directly the expression (|34|) of [io(A)] as a 
function of the rescaled variable fi in (f5U|) . Keeping terms up to the order 1/ln 2 N, one gets, using also (|57|) . 

e<m * 4_ A" ^e~» (l + ^M-lnhUV + ^-l 1 f In, - In In N + lE - ll 2 \ 
T(J3)J \ h InTV 2f InTV J ) 




The integrals of each term can be computed using Igg]) . One gets 



G(/3) _ 



1 



In' 3 AT 



/3 (Vtfi + l) 



21n 2 iV V T(/3 + 2) T(/3 + 2) 



7 + / 2 + 



(A2) 



with I = lnln A 7 " — 7s + l. Taking the logarithm of (|A2[) , one obtains G((3). By expanding in powers of /3 and comparing 
with (|2"T|) . one gets the cumulants of the position of the front. We give the velocity and diffusion constant: 



VN 



D N = 



In In N -f 

K 2 1 

3 ItlN 



lnlnAr + 1 (lnlnA 7 ) 2 



1 + — 



In A 7 



In 2 N 



2 



— In In A 7 - — 
3 6 



2C(3) 



(A3) 



Note that the first correction to the leading term can be in both cases obtained by replacing in the leading term In N 
by In N + In In A 7 : v N ~ In (In N + In In N) and D N ~ (7r 2 /3)/(ln A 7 + In In A 7 ). This is reminiscent of the observation 
in figure H] that, in the generic case, the fit was better by replacing the In A 7 by In N + 3 In In N in the theoretical 
prediction for the diffusion constant. 



2. Tree statistics 



To get subleading orders for the statistics of the tree in the exponential case, one needs to generalize the discussion 
in section 1111 Bl where we derived the leading term in the large In N expansion. The central quantity is still the 
probability r p (k) that p individuals at generation .9 + 1 have exactly k ancestors in the previous generation. But 
while at leading order it was enough to consider one coalescence at each step, one needs to take into account up to n 
simultaneous coalescences when one wishes to keep terms of arbitrary order 1/ In™ N. 

One has to assign an ancestor at generation g to each individual at generation g + 1. We start from the probability 
Wi(x) given in (|40p that the parent of an individual at position x and generation g + 1 was the i-th individual of 
generation g. In the exponential model, Wi{x) does not depend on x (see flHJ). We consider p individuals of generation 
g + 1 and we note pi the number of these individuals that are descendants of the i-th individual of generation g. The 
probability distribution of the pt is 

Proba^x, . . . , PN ) = —J*—P W? ■ --Wfr. (A4) 

Pi- • • Pn- 

One now averages over the positions of individuals at generation g, and r p (k) is simply the probability that there are 
exactly k non zero Pi's. After relabeling the individuals at generation g, one gets 

p »<*) = C0 £ sr^i^^^'"^- (A5) 

pi>i,— ,pis>i 

It is actually convenient to call n the number of Pi that are strictly larger than 1 and to write r p {k) as a sum over n: 
after another relabeling, 

^ = ( N k) E (f) E ^—^-l P , AT • • • W^W n+i .-W k ). (A6) 

Indeed, as we shall see, each term in the sum over n gives a contribution of order l/ln™ N in the final result. The 
averaged term can be expressed using the probability Wi given in (|4"Tj) : 

I" 00 r°° pPlVl-\ hPnVn+Vn+l^ VVk 

^:n Pn = (Wr---WtW n+1 ---W k )=J o dyre-yi...^ dy N e~y« {eyi+ ... + eVN)P ■ (A7) 

The technique to evaluate the integrals involved here is essentially the same as in section. Mil We first use the 
standard representation (|29|) for the denominator in the integrand. Then the integral over yi may be expressed with 
the help of the functions I p (X) defined in (|4"5|) : 



JpXn Pn = (pTT)j / dXX"- 1 I Pl (X)---IpAX)h(X) k ~ n Io(X) N - k - (A8) 



As before, for large TV, the term Iq(X) n make the integral (|A8|) dominated by values of A of order l/[TVlnTV]. It 
is sufficient to use the leading order (jUJ) for the I p (X) as next orders in A would generate terms of order 1/iV, 
which we discard throughout. Making the change of variables fj, = ATVlnTV (see (|33f ). and using the fact that 
Pi + ■ ■ ■ + p n = p — k + n, one gets for the integrand of 



A 1 Pl 1 Pn. 1 l A 



f-k „ (pi-2)!---(p n - 2)! i-Jk _^_, 
Nk-Hn^N 



1 



(k — n — /it) (In In TV — ln/i — 7^) — /i 
InTV 



can then be evaluated using 



One gets 



TPl,—,Pn 
J p,k,n 



(p 1 -2)!-..(p n -2)! (fc-1)! 



(P-1)I 



TV fe In™ TV 



r'(fc) 
r(fc) 



7B - In In TV - (fc - 1) 



In TV 



(A9) 



(A10) 



as expected, j Pl k,n has an amplitude proportional to 1/ In™ TV. To compute r p {k) for k < p to order 1/ In 2 TV, one only 
needs in (|A6|) the terms n = 1 and n = 2 (the term n = gives a contribution only for k = p): 
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fc! I (p-fc + 1)! p-' 1 ' 1 
After some algebra, one gets, for k < p, 
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(A12) 



(we used, among other things, T'(k)/T(k) + 7_e = 1 + 5 
We can now compute the (Tfc). From the recurrence 



1 

fe-l ■ 



(T p ) = l + j2r p (k){T k 



we get, using J2k r p(k) = 1 an d = °' 



fc=i 



TCfe=i r p(k) 



For the first values of p, we obtain 



(T 2 ) = lnTV + lnlnTV + o(l) 



(A13) 



(A14) 



(T 3 ) = -(InTV + lnlnTV) + o(l) 

<7 1 4> = §(lnTV + lnlnTV)-^+ (l). 



(A15) 



APPENDIX B: THE PARISI BROKEN REPLICA SYMMETRY 



The replica trick is a powerful approach to calculate the typical free energy of a sample in the theory of disordered 
systems. In the replica trick, one considers n replicas of the same random sample, one averages the product of their 
partition functions and at the end of the calculation one takes the limit n — > 0. In some cases, the n-depcndence of 
this averaged product is simple enough for the analytic continuation n — > to be unique leading to the desired free 
energy. 

In the case of mean-field spin-glasses, the situation is more complicated: the symmetry between the replicas gets 
broken as n takes non-integer values (n < 1) and remains broken in the limit n — ► 0. In this appendix we recall the 
statistical properties of the trees predicted by the Parisi theory of the broken replica symmetry [4l|, [H], [4^, [13] . 

One starts with an integer n = n$ number of replicas. These replicas are grouped into no/ni groups of ri\ replicas. 
Each of these groups of n\ replicas is decomposed into n\jni groups of 712 replicas and so on: each group of replicas 



is formed of n^/rii+i groups of rij+i replicas each. When this hierarchy consists of k levels, it is characterized by k + 1 
integers 



n = n > n\ > n 2 > ... > n k = 1. 



(Bl) 



At level i, there are a total of n/rii groups of size rii. Therefore, the probability that m distinct individuals chosen 
at random belong to the same group at level i (without specifying whether they belong or not to the same group at 
level i + 1) is 



7t~ (m) _ n(n,j - l)(n. t - 2) ■ ■ • (rtj - m + 1) 
C) ~ n(n-l)---(n-m + l) 



(B2) 



One can also associate a tree to each choice of m replicas: the m replicas are at the bottom of the tree and when two 
replicas belong to the same group at level i, but to different groups at level i + 1, their branches merge at level i. 

The various possible trees which might occur for three replicas or four replicas are shown in tables |TT] and IIIII with 
their probabilities. For example for the first tree of table HH the probability that two branches merge at level j and 
the remaining branches merge at level i is 



n(rii - n i+ i){nj - n j+ i) 
n(n - l) (n - 2) 



(B3) 



as there are n possible choices for the leftmost replica, Hi — n^+i choices for the rightmost replica and rij — nj+i choices 
for the replica at the center of the figure. The degeneracy factor is simply the number of different ways of permuting 
the roles of the replicas at the bottom of the tree. 





n(rn — m+i)(nj — n j+1 ) 
n(n - l)(n - 2) 

n(m — n i+1 )(rii — 2n i+ i) 



3 cases 



1 case 



n(n - l)(n - 2) 

TABLE II: All possible trees of three replicas, their probabilities and degeneracies. 
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1 case 



TABLE III: All possible trees of four replicas, their probabilities and degeneracies. 

In the Parisi ansatz, all the calculations are done as if all the rij's and all the ratios rij/rij+i were integers. At the 
end of the calculation, however, one takes the limit n — > and one reverses the inequality (|Bip into 



n = no < n± < 712 < ■ • • < n>k = 1- 



(B4) 



One then takes a continuous limit (k — > oo) where rii becomes a continuous variable x 

rii = x. (B5) 

In the spin glass theory [iil. l4o|. there is an ultrametric distance between pairs of replicas, related to the overlap q a ,p- 
(The distance is a decreasing function of the overlap.) This overlap q a ,p depends on the level at which the branches 
of these two replicas merge: this means that at each level i of the hierarchy, one associates a value qi of the overlap 
and that q a ^ = q t if the two replicas a and f3 belong to the same group at level i and to different groups at level i + 
(qi is an increasing function of i with qo = and q^ = 1.) In the limit k — > oo, when the rii become a continuous 
variable (|B5[) . the overlap qi becomes a increasing function q{x) = q(rii) = qi with q(0) = and q(l) = 1. 
The probability that two replicas have an overlap q ai p < qi is 

lli — fl ' 

Proba(<7„ ;/3 = q ) + Proba(g aj/3 = <7i) H h Proba(g ai/3 = qi-i) = 1 - Qafai) = -— p (B6) 

Therefore, in the n — > limit, the probability -P(q) that the overlap g ai| g between two replicas a and /3 takes the value 
q is then given by 

] P{q') dq' = lim (1 - Q 2 ) = X (B7) 

Jo 

and this leads to the famous relation (4lj | between the function q(x) and the probability distribution of the overlap 

In our models, the coalescence time between a pair of individuals in the population defines, clearly, an ultrametric 
distance. In order to see whether the statistics predicted by the replica approach remain valid for the trees of the 
exponential model discussed in the present paper, one needs to relate the overlap q(x) or the parameter x (which 
indexes the height of the hierachy) to the coalescence time T by a function T(x). It turns out that this can be achieved 
by identifying the probability e~ T that the coalescence time between two individuals is larger than T (see R 2 {T) in 
(f63l) ) with the probability that two replicas belong to different groups at level i. In other words 

e- T = l-Q a = =*P#, (B9) 
n(n — 1) 

which leads in the n — > limit to 

e- T = x. (BIO) 

With this identification, if one assumes that the statistics of the trees are given by Parisi's theory, one can compute 
all the statistical properties of the coalescence times of trees. For example, by taking the n — > limit of (|B2p , one 
gets that the probability Q m that m individuals have a coalescence time T m < T is given by 

„ r(m — x) ,„ 

^ - (Bll) 



(m-l)ir(l-a?) 
which, by taking the derivative with respect with T, gives 



I 1 dx T (x)P *9™= j°° dT TP j L T ) (B12) 



dx J dT (m- 1)! T(l - e" T ) 



This coincides with the result of the direct calculation (J62J) of the moments of the T m and shows that the statistics of 
the trees in the exponential model are the same as the ones predicted by the mean field theory of spin glasses. 



APPENDIX C: THE NEUTRAL MODEL 



In thi s appe ndix we recall some well known results on the statistical properties of the coalescence times in neutral 
models (36f E3| and derive 



We consider a population of fixed size N with non oveiiaping generations. Each individual i at a given generation 
g has ki{g) offspring at the next generation. We assume that the ki(g) are random and independent, and we call pk 
be the probability that ki{g) = k. The total number M of offspring is therefore given by 



N 



M = J2h- (Cf) 



i=l 

To keep the size of the population constant we choose N individuals at random among these M individuals. 
The probability q n that n individuals have the same parent at the previous generation is 

„ / E i fl) \ / Sjfci(fci-i)---(fci-» + i) \ (C9 , 

qn ~\ O / \ M(M-l)...(M-n + l) /■ 

For a population of large size, if pk decays fast enough with k for the moments of k to be finite, the law of large 
numbers gives that the denominator is approximatively equal to (N(k)) n and 

a - 1 / r(fc + 1) \ (C3) 
9,1 " iV™-i(fc)« \r(fe-n + l)/" [ ' 

We see (when the moments of k are finite) that qi is much larger than all the other q n when the size N of the 
population is large, and therefore in the ancestry of a finite number n of individuals, branches coalesce only by pairs. 
Similarly, the probability that two or more pairs of individuals coalesce at the same generation is negligible. 

Let T n (g) be the age of the most recent common ancestor of a group of n individuals at generation g. As for large N 
only coalesences by pairs may occur from one generation to the previous one, one has 

{T„(q) + 1 with probability 1 — in(n — l)qo 
, (C4) 
T„_i(g) + 1 with probability ±n(n — 1)<?2- 

In the steady state [HI, this implies that 

= (l - ^^4) <(1 + T n y) + ^j^q 2 ((l + T n . x f) (C5) 
and using the fact that T\{g) = and one gets 

(T„)=(2-£)I. (06) 

We see that all the times T n scale like N (since qi ~ iV -1 ) and that 

(T 3 ) = 4 (T 4 ) = 3 (T n ) = 2(n - f ) 

(T 2 ) 3' (T 2 ) 2' (T a ) n ' 

One can also calculate from (|C5p higher moments of the T n 's or their generating functions 

m?) n ((r 3 ) 2 ) is 



(T 2 )2 ' (T £ 



(07) 



(08) 



These distributions of the T„ as well as their correlations are universal (in the sense that they do not depend on the 
details of the distribution of the pk 's) . 



[1] L. Peliti, Lectures at the Summer College on Frustrated System, Trieste (1997), cond-mat/9712027. 
[2] D. A. Kessler, H. Levine, D. Ridgway, and L. Tsimring, J. Stat. Phys. 87, 519 (1997). 
[3] L. Tsimring, H. Levine, and D. A. Kessler, Phys. Rev. Lett. 76, 4440 (1996). 
[4] P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993). 
[5] M. Kloster and C. Tang, Phys. Rev. Lett. 92, 038101 (2004). 



[6] M. Kloster, Phys. Rev. Lett. 95, 168701 (2005). 
[7] R. E. Snyder, Ecol. 84, 1333 (2003). 

[8] E. Bmnet, B. Derrida, A. H. Mueller, and S. Munier, Phys. Rev. E 73, 056126 (2006). 
[9] E. Brunet, B. Derrida, A. H. Mueller, and S. Munier, Europhys. Lett. 76, 1 (2006). 
[10] R. A. Fisher, Annals of Eugenics 7, 355 (1937). 

[11] A. Kolmogorov, I. Petrovsky, and N. Piscounov, Bull. Univ. Etat Moscou, A 1, 1 (1937). 

[12] W. van Saarloos, Phys. Rep. 386, 29 (2003). 

[13] B. Derrida and H. Spohn, J. Stat. Phys. 51, 817 (1988). 

[14] E. Bmnet and B. Derrida, Phys. Rev. E 70, 016106 (2004). 

[15] H.-P. Breuer, W. Huber, and F. Petruccione, Europhys. Lett. 30, 69 (1995). 

[16] C. R. Doering, C. Mueller, and P. Smereka, Physica A 325, 243 (2003). 

[17] A. Lemarchand and B. Nowakowski, Journal of Chemical Physics 111, 6190 (1999). 

[18] E. Moro, Phys. Rev. Lett. 87, 238303 (2001). 

[19] P. L. Krapivsky and S. N. Majumdar, Phys. Rev. Lett. 85, 5492 (2000). 
[20] E. Iancu, A. H. Mueller, and S. Munier, Phys. Lett. B 606, 342 (2005). 
[21] S. Munier and R. Peschanski, Phys. Rev. Lett. 91, 232001 (2003). 
[22] A. H. Mueller and A. I. Shoshi, Nucl. Phys. B 692, 175 (2004). 
[23] C. Mueller and R. B. Sowers, J. Funct. Anal. 128, 439 (1995). 
[24] D. Panja, Phys. Rep. 393, 87 (2004). 

[25] J. G. Conlon and C. R. Doering, J. Stat. Phys. 120, 421 (2005). 
[26] C. Escudero, Phys. Rev. E 70, 041102 (2004). 
[27] E. Moro, Phys. Rev. E 69, 060101(R) (2004). 

[28] J. Mai, I. M. Sokolov, and A. Blumen, Phys. Rev. Lett. 77, 4462 (1996). 

[29] D. A. Kessler, Z. Ner, and L. M. Sander, Phys. Rev. E 58, 107 (1998). 

[30] E. Moro, Phys. Rev. E 70, 045102(R) (2004). 

[31] L. Pechenik and H. Levine, Phys. Rev. E 59, 3893 (1999). 

[32] A. Lemarchand, A. Lesne, and M. Mareschal, Phys. Rev. E 51, 4457 (1995). 

[33] H. P. McKean, Comm. Pure Appl. Math. 28, 323 (1975). 

[34] M. D. Bramson, Communications In Pure and Applied Mathematics 31, 531 (1978). 

[35] J. F. C. Kingman, Stoch. Proc. Appl. 13, 235 (1982). 

[36] J. F. C. Kingman, J. Appl. Probab. 19A, 27 (1982). 

[37] J. Pitman, Ann. Probab. 27, 1870 (1999). 

[38] J. Schweinsberg, Elect. Journ. Prob. 5, 1 (2000). 

[39] E. Bolthausen and A.-S. Sznitman, Com. Math. Phys. 197, 247 (1998). 

[40] S. Tavare, D. J. Balding, R. C. Griffiths, and P. Donelly, Genetics 145, 505 (1997). 

[41] G. Parisi, Phys. Rev. Lett. 50, 1946 (1983). 

[42] E. Brunet and B. Derrida, Phys. Rev. E 56, 2597 (1997). 

[43] E. Brunet and B. Derrida, J. Stat. Phys. 103, 269 (2001). 

[44] M. Mezard, G. Parisi, N. Sourlas, G. Toulouse, and M. A. Virasoro, Journal de Physique 45, 843 (1984). 

[45] M. Mezard, G. Parisi, and M. A. Virasoro, Spin glass theory and beyond (World Scientific Lecture Notes in Physic, 1987). 

[46] D. Ruelle, Com. Math. Phys 108, 225 (1987). 

[47] G. Parisi, J. Phys. A 13, 1101 (1980). 

[48] D. Simon and B. Derrida, J. Stat. Mech., P05002 (2006). 



